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Large-eddy  simulation  of  conductive  flows  at  low 
magnetic  Reynolds  number 

By  B.  Knaepen  and  P.  Moin 


1.  Introduction 

In  large-eddy  simulations  (LES),  only  the  large-scale  structures  of  the  flow  are  simu- 
lated directly  while  the  small-scale  structures  are  taken  into  account  through  a model, 
referred  to  as  the  subgrid-scale  (SGS)  model.  This  results  in  significant  reduction  in  com- 
putational requirements.  The  trade-off  comes  in  the  extra  modelling  effort  that  has  to  be 
produced  in  order  to  adequately  take  into  account  the  discarded  small  scales  structures. 
In  the  computational  fluid  mechanics  community,  the  most  widely  used  SGS  model  is 
the  Smagorinsky  eddy  viscosity  model  (Smagorinsky  1963).  This  model  has  gained  an 
even  bigger  practical  interest  following  the  work  of  Germano  et  al.  (1991)  in  which  the 
dynamic  procedure  was  introduced.  The  dynamic  procedure  enables  optimization  “on  the 
fly”  of  the  arbitrary  scaling  factor  that  is  inherently  present  in  the  original  Smagorinsky 
model  (more  details  below)  and  thus  allows  the  model  to  automatically  adapt  to  the  flow 
being  studied. 

In  this  paper  we  study  the  LES  method  with  dynamic  procedure  in  the  context  of 
conductive  flows  subject  to  an  applied  external  magnetic  field  at  low  magnetic  Reynolds 
number  Rm-  These  kind  of  flows  are  encountered  in  many  industrial  applications.  For 
example,  in  the  steel  industry,  applied  magnetic  fields  can  be  used  to  damp  turbulence  in 
the  casting  process.  In  nuclear  fusion  devices  (Tokamaks),  liquid-lithium  flows  are  used  as 
coolant  blankets  and  interact  with  the  surrounding  magnetic  field  that  drives  and  confines 
the  fusion  plasma.  Also,  in  experimental  facilities  investigating  the  dynamo  effect,  the 
flow  consists  of  liquid-sodium  for  which  the  Prandtl  number  and,  as  a consequence,  the 
magnetic  Reynolds  number  is  low. 

Most  of  the  previous  works  considering  LES  in  the  case  of  MHD  flows  have  been 
directed  towards  flows  at  high  Rm,  where  the  full  non-lineax  MHD  equations  have  to  be 
used,  or  towards  flows  in  the  absence  of  an  external  magnetic  field:  Yoshizawa  (1987); 
Zhou  & Vahala  (1991);  Theobald  et  al.  (1994);  AguUo  et  al.  (2001);  Muller  & Carati 
(2002).  To  our  knowledge,  the  only  attempt  to  study  MHD  turbulence  at  low  magnetic 
Reynolds  number  from  the  LES  point  of  view  is  due  to  Shimomura  (1991).  In  that 
work,  the  author  considers  the  case  of  magnetohydrodynamic  turbulent  channel  flow  and 
introduces  a new  SGS  model  designed  to  incorporate  the  effects  of  the  applied  uniform 
magnetic  field.  However,  like  the  original  Smagorinsky  model,  the  model  of  Shimomura 
contains  constants  that  have  to  be  adjusted  and  are  thus  flow  dependent.  In  the  present 
paper,  we  show  that  it  is  possible  to  circumvent  this  problem  by  the  use  of  the  dynamic 
Smagorinsky  model. 

Our  attention  is  focused  here  on  the  case  of  homogeneous  (initially  isotropic)  de- 
caying turbulence.  The  numerical  simulations  performed  mimic  the  thought  experiment 
described  in  Moffatt  (1967)  in  which  an  initially  homogeneous  isotropic  conductive  flow 
is  suddenly  subjected  to  an  applied  magnetic  field  and  freely  decays  without  any  forcing. 
Note  that  this  flow  was  first  studied  numerically  by  Schumann  (1976).  It  is  well  known 
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that  in  that  case,  extra  damping  of  turbulence  occurs  due  to  the  Joule  effect  and  that 
the  flow  tends  to  become  progressively  independent  of  the  coordinate  along  the  direction 
of  the  magnetic  field.  Our  comparison  of  filtered  direct  numerical  simulation  (DNS)  pre- 
dictions and  LES  predictions  show  that  the  dynamic  Smagorinsky  model  enables  one  to 
capture  successfully  the  flow  with  LES,  and  that  it  automatically  incorporates  the  effect 
of  the  magnetic  field  on  the  turbulence. 

Our  paper  is  organized  as  follows.  In  the  next  section  we  summarize  the  LES  approach 
in  the  case  of  MHD  turbulence  at  low  Rm  and  recall  the  definition  of  the  dynamic 
Smagorinsky  model.  In  Sec.  3 we  describe  the  parameters  of  the  numerical  experiments 
performed  and  the  code  used.  Section  4 is  devoted  to  the  comparison  of  filtered  DNS 
results  and  LES  results.  Conclusions  are  presented  in  Sec.  5. 


2.  LES  at  low  magnetic  Reynolds  number 

For  homogeneous  flows,  the  magnetic  Reynolds  number  can  be  defined  through  the 
following  relation: 


Rm  - 


uL 

V 


(2.1) 


where  u = (uiUi)  /3  is  the  r.m.s.  of  the  fluctuating  velocity  Ui,  L is  the  integral  length 
scale  of  the  flow  and  rj  is  the  magnetic  diffusivity.  Rm  represents  the  relative  importance 
of  the  non-linear  terms  and  the  diffusion  term  in  the  magnetic  induction  equation.  In  the 
limit  of  low  Rm,  the  MHD  equations  can  be  simplified  considerably  Roberts  (1967).  It  is 
indeed  possible  to  close  independently  the  momentum  equation  and  take  into  account  the 
effect  of  the  magnetic  field  through  an  extra  damping  term.  This  approximation  is  known 
as  the  quasi-static  ( QS)  approximation  and  along  with  the  incompressibility  condition, 
diUi  = 0,  reads. 


dtUi  — —di{p/p)  - UjdjUi  — 


■n 


A -t- 


(2.2) 


where  p is  the  sum  of  the  kinematic  and  magnetic  pressures,  p is  the  fluid  density,  i/  the 
kinematic  viscosity,  is  the  applied  external  magnetic  field  and  A~^  is  the  inverse 
of  the  Laplacian  operator.  Note  that  hris  by  convention  been  eiligned  with  the  z- 
direction  and  expressed  in  Alfven  units. 

LES  equations  are  obtained  by  filtering  (2.2).  The  filtered  velocity  Ui  is  defined  by. 


Wi(x)  = /(7(x,y)ui(y)dy,  (2.3) 

where  G is  a smoothing  kernel  that  eliminates  the  small  scale  part  of  Ui  and  satisfies  the 
relation. 


J G(x,y)dy  = 1. 

In  terms  of  the  filtered  velocity,  Eq.  2.2  can  be  written  as. 


(2.4) 


dtUi  = -di{p/p)  — UjdjUi  - 


A ^d;:d^Ui  + i/Aui-djTij, 


V 


(2.5) 


where, 
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T ij  — UjXLj 


(2.6) 


is  the  subgrid-scale  stress  tensor.  In  order  to  close  (2.5),  Tij  has  to  be  expressed  only  in 
terms  of  the  filtered  velocity. 

Is  is  interesting  to  note  that  does  not  depend  explicitly  on  the  magnetic  field. 
Indeed,  the  magnetic  contribution,  being  a linear  term  in  (2.2),  does  not  require  an 
explicit  SGS  counterpart  in  the  LES  equation.  In  the  case  of  non-conductive  flows,  the 
most  widely  used  model  for  ry  is  the  Smagorinsky  model: 

Tij  = -2CsA'^\S\Sij,  Sij  = ^(diUj  + djUi),  S = yj2SijSij,  (2.7) 


where  A is  the  width  of  the  filter  G and  Ca  is  the  Smagorinsky  constant.  However,  as 
was  observed  by  Shimonura  in  Shimomura  (1991),  model  (2.7)  with  Cg  optimized  for  a 
non-conductive  flow,  is  not  adequate  in  the  present  context.  This  is  easily  understood 
when  one  refers  to  the  fact  that  the  applied  external  magnetic  field  has  the  tendency 
to  suppress  non-linear  transfers  in  the  velocity  field.  Thus,  the  effect  of  the  SGS  stress 
tensor  needs  to  be  decreased  for  conductive  flows  when  the  external  magnetic  field  is 
switched  on.  This  should  result  in  a lower  optimal  value  for  the  Smagorinsky  constant. 

Contrary  to  what  was  done  in  Shimomura  (1991),  we  will  not  explicitly  modify  the 
Smagorinsky  model  to  incorporate  the  eflFect  of  the  magnetic  field  but  rather  we  will 
make  use  of  the  d3mamic  procedure  (Germano  et  al.  1991)  to  optimize  the  value  of  Ca- 

To  define  the  dynamic  procedure  one  introduces  a second,  coarser  filter  called  the 
test-filter  which  we  denote  by  TT..  Applying  this  filter  in  addition  to  (2.5)  yields. 


dtUi  = -di{p/p)  - UjdjUi  - S.  . A'^d^dgUi  + vAui  - djfij  - djLij, 


(2.8) 


where  Lij  = UiUj  — UiUj  is  the  Leonard  tensor  which  does  not  require  any  modeling  since 
it  is  expressed  in  closed  form  using  the  filtered  velocity  u,.  The  sum  Tij  -t-  Ljj  represents 

the  SGS  stress  tensor  Tjj  of  the  combined  ~ filter  and  we  have  the  well-known  Germano 
identity, 

Lij  = Tij -Tij.  (2.9) 

Assuming  that  Ui  and  Uj  are  self-similar,  a suitable  model  for  Tij  should  be, 

fij  = -2Cat^(E\%j,  %j  = ^{di%  + dj%i),  (2.10) 

where  A is  the  width  of  the  ~ filter.  When  one  models  Tij  and  Tij  using  (2.7)  and 
(2.10),  the  Germano  identity  will  undoubtedly  be  violated.  However,  the  constant  Ca 
can  be  chosen  in  such  a way  as  to  minimize  the  difference  (in  the  least  square  sense) 
between  both  sides  of  (2.9).  Assuming  homogeneity  in  all  directions,  the  optimal  choice 
is  (Lilly  1992): 
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where  (•  • • ) denotes  spatial  averaging.  Thus,  although  model  (2.7)  does  not  explicitly 
incorporate  the  magnetic  field,  the  dynamic  constant  Cg  should  adjust  to  an  appropriate 
value  if  the  scale-similarity  hypothesis  is  at  all  justified. 


3.  Numerical  experiments 

To  assess  the  dynamic  Smagorinsky  model  in  the  present  context,  we  have  built  a 
set  of  DNS  databases.  Since  we  restrict  our  attention  to  a cubic  domain  with  periodic 
boundary  conditions,  a spectral  (dealiased)  code  has  been  used. 

The  velocity  field  can  then  be  initialized  in  Fourier  space.  The  initial  mode  amplitudes 
are  computed  to  match  the  spectra  of  the  Comte-Bellot  & Corrsin  (1971)  experiment  at 
stage  1 (see  Rogallo  (1981)  for  more  information  on  the  procedure).  Phases  are  initially 
random  and  the  flow  is  left  to  freely  decay  until  the  skewness  of  the  velocity  derivative 
reaches  a quasi  constant  value  of  5 = —0.4.  At  this  time  (hereafter  referred  to  as  to) 
the  flow  is  considered  “physical”  and  used  as  the  initial  condition  for  all  of  our  runs.  All 
the  DNS  simulations  are  done  in  a {2ir)^  box  using  a resolution  of  512^  Fourier  modes 
and  the  viscosity  is  set  to  v = 0.006.  Other  relevant  quantities  measured  at  t = to  are 
summarized  here: 


uL 

Re  = — = 380  (integral  Reynolds  number), 

/ j5 

R\  = y ~ (microscale  Reynolds  number), 
3 V? 

Teddy  = X — = 0.238  (eddy  turnover  time). 

^ € 


(3.1) 

(3.2) 

(3.3) 


Three  test  cases  have  been  considered.  The  first  one  corresponds  to  a decaying  flow 
without  the  addition  of  any  applied  external  magnetic  field.  The  other  two  numerical 
experiments  are  distinguished  by  their  characteristic  interaction  numbers  at  to-  The  in- 
teraction number  N (also  known  as  the  Stuart  number)  is  defined  as  follows: 


N = 


(Br^f  L 

77  u 


(3.4) 


N measures  the  relative  strengths  of  the  magnetic  damping  term  and  the  non-linear  term 
in  (2.2).  Two  cases  are  examined  here,  iV  = 1 and  N = 10. 

To  illustrate  the  effect  of  the  external  magnetic  field,  the  time  history  of  the  kinetic 
energy  density  E = ^ f dx^Ui(x)ui(x)  of  the  flow  for  the  three  cases  is  presented  in  Fig. 
1.  As  in  all  the  following  figures,  time  has  been  normalized  by  the  initial  eddy  turnover 
time  and  the  non-dimensional  time  is  denoted  by  t*. 


4.  LES  results 

In  order  to  assess  the  LES  method,  some  DNS  snapshots  of  the  flow  field  have  been 
truncated  from  the  512^  resolution  to  a 32^  resolution  (using  a sharp  Fourier  cut-off). 
The  initial  condition  for  the  LES  runs  has  been  obtained  in  the  same  way  by  truncating 
the  DNS  field  at  t = to-  As  can  be  easily  seen  from  (2.11),  the  only  free  parameter  in  the 

subgrid-scale  model  is  the  ratio  A/ A.  As  is  standard  practice  in  LES  of  non-conductive 
flows,  the  value  of  2 is  adopted  here. 
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Figure  1.  Time  history  of  the  kinetic  energy  density  E.  The  solid  line  represents  the  flow 
decay  without  applied  magnetic  field,  the  deished  curve  corresponds  to  the  case  JV  = 1 and  the 
dash-dot  curve  corresponds  to  the  case  N = 10.  The  dotted  line  represents  the  time  to  at  which 
the  magnetic  field  is  switched  on. 


Figure  2.  Time  history  of  the  resolved  kinetic  energy  density:  LES  vs.  filtered  DNS.  The  solid 
lines  represent  the  predictions  of  the  LES.  The  diamonds  represent  the  filtered  kinetic  energy 
density  obtained  by  truncating  the  DNS  fields  to  a 32®  resolution.  The  dotted  lines  represent 
the  predictions  of  the  unresolved  DNS  (performed  on  the  32®  mesh)  without  any  subgrid-scale 
model. 


4.1.  Kinetic  Energy 

Figure  2 represents  the  time  evolution  of  the  resolved  kinetic  energy  density  Er  = 
^ f dx^Ui(x)ui(x)  predicted  by  the  LES  and  compared  to  the  filtered  DNS.  On  each 
plot,  a third  curve  representing  a simulation  (referred  to  as  ‘unresolved  DNS’)  on  the 
32^  mesh  without  SGS  modeling  is  added  to  stress  the  action  of  the  subgrid-scale  model. 
The  case  5®®*  = 0 serves  as  a benchmark  to  check  that  in  the  case  of  non-conductive 
flows,  our  LES  code  behaves  as  expected.  In  both  the  cases  N = 1 and  N = 10,  the  LES 
performs  remarkably  well.  In  the  case  N = 1 the  difference  between  LES  and  unresolved 
DNS  is  very  clear.  In  the  case  N = 10  and  for  this  diagnostic,  the  unresolved  DNS  does 
not  depart  significantly  from  the  filtered  DNS  and  LES. 

4.2.  Energy  spectra 

The  energy  spectrum  constitutes  a finer  diagnostic  than  the  kinetic  energy  density  since 
it  retains  information  about  the  repartition  of  energy  among  the  different  scales  of  the 
flow.  Figures  3 and  4 contain  respectively  kinetic  energy  spectra  for  the  cases  N = 1 
and  N = 10  at  several  instants  in  the  simulation.  Agreement  between  LES  and  DNS  is 
again  excellent.  The  unresolved  DNS  exhibits  the  usual  pile  up  of  energy  at  high  wave 
numbers  that  results  from  lack  of  resolution.  Thus,  although  in  terms  of  global  energy 
the  unresolved  DNS  does  perform  quite  well  in  the  IV  = 10  case,  the  spectral  properties 
of  the  flow  are  poorly  predicted  and  the  LES  does  a much  better  job. 
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Figure  3.  Three-dimensional  kinetic  energy  spectra  for  the  case  N — 1.  The  solid  line  represents 
the  filtered  DNS,  the  dashed  line  represents  the  LES  and  the  dotted  line  corresponds  to  the 
unresolved  DNS.  The  times  at  which  the  spectra  are  calculated  are  indicated  in  the  plots. 


k k k 


Figure  4.  Three-dimensional  kinetic  energy  spectra  for  the  case  N = 10.  The  solid  line  repre- 
sents the  filtered  DNS,  the  dashed  line  represents  the  LES  and  the  dotted  line  corresponds  to 
the  unresolved  DNS.  The  times  at  which  the  spectra  are  calculated  are  indicated  in  the  plots. 


4.3.  Dissipation  rate  budget 

The  evolution  of  the  resolved  kinetic  energy  density  Er  is  governed  by, 

dEn  1 f 

~dT  J (4.1) 

where, 

/’Dext\2  

€i/  “ ^l^SijSij^  ~ — U{^  ^sgs  ~ rijS^j,  (4.2) 

Prom  Fig.  2 it  is  clear  that  the  LES  predicts  the  total  dissipation  rate  very  well.  It  is 
however  instructive  to  know  how  the  different  contributions  (4.2)  are  reproduced  sepa- 
rately. 

In  Fig.  5 we  present  the  time  evolution  of  the  resolved  viscous  dissipation  rate  for  the 
LES,  the  filtered  DNS  and  the  unresolved  DNS.  The  figures  show  that  the  LES  reproduces 
this  diagnostic  very  well,  whereas  the  unresolved  DNS  systematically  overestimates  it. 
This  results  from  the  pile-up  of  energy  near  the  LES  cut-off  where  most  of  the  resolved 
viscous  dissipation  occurs. 

Figure  6 represents  the  time  evolution  of  the  resolved  magnetic  dissipation  eg-  Again 
the  LES  predictions  match  the  filtered  DNS  results  very  well.  The  differences  with  the 
unresolved  DNS  are  less  severe  than  for  the  viscous  dissipation.  This  is  to  be  expected 
since  the  magnetic  dissipation  occurs  at  every  scale  in  the  flow  and  its  overall  intensity 
is  thus  less  contaminated  by  the  pile-up  of  energy  occurring  near  the  LES  cut-off  for  the 
unresolved  DNS.  It  is  also  interesting  to  note  that  in  the  case  N = 10,  the  magnetic 
dissipation  falls  very  rapidly  with  time  at  the  early  stages  of  the  decay.  This  happens 
because  in  that  case  the  flow  quickly  becomes  fairly  independent  of  the  z-direction  which 
is  parallel  to  the  magnetic  field  (this  is  illustrated  further  in  Section  4.4). 

Finally,  in  Fig.  7 the  evolution  of  the  subgrid-scale  transfer  rate  Csga  is  presented.  For 
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Figure  5.  Resolved  viscous  dissipation  rate  for  the  cases  N = I (left)  and  N = 10  (right). 
In  each  plot,  the  solid  line  represents  the  LES,  the  diamonds  represent  the  filtered  DNS  and  the 
dotted  line  represents  the  unresolved  DNS. 


Figure  6.  Resolved  magnetic  dissipation  rate  es  for  the  cases  N = 1 (left)  and  N = 10  (right). 
In  each  plot,  the  solid  line  represents  the  LES,  the  diamonds  represent  the  filtered  DNS  and  the 
dotted  line  represents  the  unresolved  DNS. 

the  LES,  6ggs  is  obtained  by  substituting  the  dynamic  Smagorsinky  in  place  of  Ty  in 
(4.2).  DNS  results  are  obtained  by  computing  the  exact  = vfuf  — UiUj  from  the  DNS 

data  fields.  Obviously,  no  (non-zero)  results  are  available  for  the  unresolved  DNS.  Prom 
the  plots  we  see  that  the  subgrid-scale  transfer  rate  is  initially  overestimated  in  the  LES. 
This  is  not  a surprise  since  the  dynamic  procedure  usually  needs  a httle  time  to  settle. 
Soon  after  this  short  transient  time,  the  agreement  between  LES  and  DNS  predictions 
is  very  good. 


4.4.  Flow  structures 

It  is  well  known  that  the  extra  damping  term  present  in  the  quasi-static  approximation 
(2.2)  leads  to  a progressive  suppression  of  spatial  variations  in  the  flow  along  the  direction 
of  the  magnetic  field.  It  is  thus  important  to  assess  how  well  the  LES  is  able  to  reproduce 
this  feature.  To  that  end,  we  have  plotted  in  Figures  8 and  9 the  contours  of  the  kinetic 
energy  density  at  three  different  times  respectively  for  the  filtered  DNS  and  the  LES 
(only  the  case  N = 10  is  shown  because  the  effect  is  more  pronounced  for  strong  applied 
magnetic  fields).  Since  the  turbulence  is  decaying  with  time,  the  colormap  had  to  be 
rescaled  for  each  of  the  single  plots.  However,  the  same  colormap  for  the  filtered  DNS 
and  LES  has  been  used  in  the  plots  corresponding  to  the  same  times.  As  is  obvious  from 
the  plots,  the  LES  reproduces  the  filtered  DNS  structures  very  well.  The  unresolved  DNS 
also  captures  the  large-scales  structures  of  the  flow  (no  plots  shown)  but  due  to  the  pile 


t*  = 0.285  r = 1.02  r = 2.68 

Figure  9.  Contours  of  the  kinetic  energy  obtained  from  the  LES  (AT  = 10  case).  The  different 
times  at  which  the  contours  are  calculated  are  indicated  under  the  plots  and  are  the  same  ones 
as  those  used  in  Fig.  8. 


up  of  energy  it  predicts  a significant  amount  of  small  scales;  the  corresponding  contour 
plots  contain  some  excessive  ‘noise’  on  top  of  the  large-scales  structures. 

4.5.  Dynamic  Constant 

As  recalled  in  Section  2,  the  dynamic  procedure  is  designed  to  optimize  the  scaling 
constant  Cg  present  in  the  Smagorinsky  subgrid-scale  model.  We  also  mentioned  that 
when  an  external  magnetic  field  is  present,  the  value  of  C,  should  decrease  since  non- 
linear transfers  are  reduced  and  the  effect  of  the  subgrid-scale  model  should  be  damped 
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Figure  10.  Time  history  of  the  dynsmiic  Smagorinsky  constant.  The  solid  line  represents  the 
case  with  no  magnetic  field,  the  dashed  line  corresponds  to  the  case  = 1 and  the  case  iV  = 10 
is  represented  by  the  dotted  line. 
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accordingly.  In  Fig.  10  the  evolution  of  is  plotted.  Prom  the  figure,  it  is  clear  that 
the  expected  behavior  is  observed.  As  the  interaction  number  is  increased,  the  dynamic 
procedure  automatically  adapts  the  value  of  the  Smagorinsky  constant  to  decrease  the 
effect  of  the  subgrid-scale  model. 


5.  Conclusions 

In  this  article  we  have  shown  that  the  dynamic  Smagorinsky  model  can  be  used  to 
perform  large-eddy  simulations  of  flows  subject  to  an  applied  external  magnetic  field  at 
low  magnetic  Reynolds  number,  Rm-  Although  this  subgrid-scale  model  was  not  designed 
for  this  appHcation,  its  behavior  is  excellent  owing  to  its  adaptation  to  the  flow  and  the 
applied  magnetic  field  through  the  dynamic  procedure.  The  model  can  be  considered 
robust  since  it  works  equally  well  for  interaction  (or  Stuart)  numbers  ranging  firom  N = 0 
(no  magnetic  field)  to  iV  = 10  (for  which  the  flow  becomes  nearly  two  dimensional).  In 
the  future,  the  same  model  will  be  tested  in  more  complex  geometries. 
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